﻿#include "DataSourceImageGDAL.h"
#include <gdal_priv.h>
#include "geoimage.h"
#include <qDebug>
#include <ogr_spatialref.h>
#include <math.h>
#include "CHuGlobalMercator.h"


using namespace rxtiler;

DataSourceImageGDAL::DataSourceImageGDAL()
{

}

DataSourceImageGDAL::~DataSourceImageGDAL()
{
    this->close();
}

bool DataSourceImageGDAL::close()
{
    if(poDataset){
        GDALClose(poDataset );
        poDataset = 0;
    }
    return true;
}
bool DataSourceImageGDAL::open()
{
    //3857
    poDataset = ( GDALDataset* )GDALOpen( _path.toStdString().c_str(), GA_ReadOnly );
    if(!poDataset) return false;
}

rxtiler::Extent DataSourceImageGDAL::getExtent()
{
    double adfGeoTransform[6];
    if (poDataset && poDataset->GetGeoTransform(adfGeoTransform) == CE_None){
        int x=poDataset->GetRasterXSize();
        int y=poDataset->GetRasterYSize();
        double left=adfGeoTransform[0];
        double top=adfGeoTransform[3];
        double right=adfGeoTransform[0]+x*adfGeoTransform[1]+y*adfGeoTransform[2];
        double bottom=adfGeoTransform[3]+x*adfGeoTransform[4]+y*adfGeoTransform[5];
        _tif.left = left;
        _tif.right = right;
        _tif.bottom = bottom;
        _tif.top = top;
        return rxtiler::Extent(left,right,bottom,top);
    }else
    {
        return rxtiler::Extent(0,0,0,0);
    }
}


osg::Image *DataSourceImageGDAL::getImage(rxtiler::TileKey &key, rxtiler::Extent &extent)
{
    if(!poDataset){
        return 0;
    }

    int img_size = this->getTileSize();
    if(!extent.cover(_tif)){
        osg::Image* img = new osg::Image();
        img->allocateImage(img_size,img_size,4,GL_RGBA,GL_UNSIGNED_BYTE);
        return img;
    }
    bool needResample= true;
    if(_tif.contain(extent)){
        needResample = false;
    }

    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();


    Extent ext2 = extent.intersect(_tif);
    int x = iWidth * (ext2.left - _tif.left) / _tif.width();
    int y = iHeight * ( _tif.top - ext2.top) / _tif.height();
    int XSize = ::ceil(iWidth * ext2.width() /  _tif.width());
    int YSize = ::ceil(iHeight * ext2.height()  /  _tif.height());
    if(x + XSize >iWidth ) XSize = iWidth - x;
    if(y + YSize >iHeight ) YSize = iHeight - y;

    int destX = this->getTileSize();
    int destY = destX;
    osg::Image* img = this->__read(x,y,XSize,YSize,destX,destY);
    if(!img) return 0;
    if(needResample){
        osg::Image* img2 = new osg::Image();
        img2->allocateImage(img_size,img_size,4,GL_RGBA,GL_UNSIGNED_BYTE);
        mapImage(img,ext2,img2,extent);
        return img2;
    }else{
        return img;
    }
}

osg::Image *DataSourceImageGDAL::getImages(TileKey &key, Extent &extent)
{
    if(!poDataset){
        return 0;
    }

    //1赋值extent
    double adfDstGeoTransform[6]={0,0,0,0,0,0};
    poDataset->GetGeoTransform(adfDstGeoTransform);
    int nXsize=poDataset->GetRasterXSize();
    int nYsize=poDataset->GetRasterYSize();



    double dbX[4] ={0,0,0,0};

    double dbY[4] = {0,0,0,0};

    dbX[0] = adfDstGeoTransform[0];    //左上角点坐标

    dbY[0] = adfDstGeoTransform[3];


    //右上角坐标

    dbX[1] = adfDstGeoTransform[0] + nXsize*adfDstGeoTransform[1];

    dbY[1] = adfDstGeoTransform[3];


    //右下角点坐标

    dbX[2] = adfDstGeoTransform[0] + nXsize*adfDstGeoTransform[1] + nYsize*adfDstGeoTransform[2];

    dbY[2] = adfDstGeoTransform[3] + nXsize*adfDstGeoTransform[4] + nYsize*adfDstGeoTransform[5];


    //左下角坐标

    dbX[3] = adfDstGeoTransform[0];

    dbY[3] = adfDstGeoTransform[3] + nYsize*adfDstGeoTransform[5];





    if(dbY[0]>85.05112878)
    {
        dbY[0]=85.05112878;
    }
    if(dbY[1]>85.05112878)
    {
        dbY[1]=85.05112878;
    }
    if(dbY[2]<-85.05112878)
    {
        dbY[2]=-85.05112878;
    }
    if(dbY[3]<-85.05112878)
    {
        dbY[3]=-85.05112878;
    }

    if(dbX[0]<-180)
    {
        dbX[0]=-180;
    }
    if(dbX[3]<-180)
    {
        dbX[3]=-180;
    }

    //左上角
    double topleft[2] = {0,0};
    topleft[0] = dbX[0];
    topleft[1] = dbY[0];

    //右下角
    double loweright [2] = {0,0};

    loweright[0] = dbX[2];
    loweright[1] = dbY[2];

    CHuGlobalMercator m_self(this->getTileSize());
    double LatLonbound5[4] ={0,0,0,0};



    int mx = (pow(2,(double)key.lod) - 1) - key.x;
//    key.x = mx;
    m_self.TileLatLonBounds(key.y,mx,key.lod,LatLonbound5);


    extent.left  =LatLonbound5[1];
    extent.top  =LatLonbound5[2];
    extent.right  =LatLonbound5[3];
    extent.bottom  =LatLonbound5[0];
    _tif.left  =adfDstGeoTransform[0];
    _tif.top  =adfDstGeoTransform[3];
    _tif.right  =adfDstGeoTransform[0] + nXsize*adfDstGeoTransform[1] + nYsize*adfDstGeoTransform[2];
    _tif.bottom  =adfDstGeoTransform[3] + nXsize*adfDstGeoTransform[4] + nYsize*adfDstGeoTransform[5];



    //3计算图片
    int img_size = this->getTileSize();
    if(!extent.cover(_tif)){
        osg::Image* img = new osg::Image();
        img->allocateImage(img_size,img_size,4,GL_RGBA,GL_UNSIGNED_BYTE);
        return img;
    }
    bool needResample= true;
    if(_tif.contain(extent)){
        needResample = false;
    }
    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();


    Extent ext2 = extent.intersect(_tif);
    int x = iWidth * (ext2.left - _tif.left) / _tif.width();

    int y = iHeight * ( _tif.top - ext2.top) / _tif.height();

    int XSize = ::ceil(iWidth * ext2.width() /  _tif.width());
    int YSize = ::ceil(iHeight * ext2.height()  /  _tif.height());

    if(x + XSize >iWidth )
    {
        XSize = iWidth - x;
    }
    if(y + YSize >iHeight )
    {
        YSize = iHeight - y;
    }

    int destX = this->getTileSize();
    int destY = destX;

    osg::Image* img = this->__reads(x,y,XSize,YSize,XSize,YSize);
    osg::Image* outimg = new osg::Image();
    outimg->allocateImage(destX,destY,4,GL_RGBA,GL_UNSIGNED_BYTE);


    double latlan1[2] ={0,0};
    double latlan2[2] ={0,0};
    for (int i = 0; i < destY; ++i)
    {
        for (int j = 0; j <destX ; ++j)
        {
            m_self.LatLonToMeters(extent.left, extent.top , latlan1);
            m_self.LatLonToMeters(extent.right, extent.bottom , latlan2);
            double m_res=  m_self.Resolution(key.lod);
            double m_lon = latlan1[0] + m_res*i;
            double m_lat = latlan1[1] -m_res*j;


            m_self.MetersToLatLon(m_lon, m_lat , latlan1);
            int x3 = XSize * (latlan1[1]  - extent.left) /  extent.width();
            int y3 = YSize * (ext2.top - latlan1[0]) / ext2.height();;

            osg::Vec4 c = img->getColor(x3,y3);
            outimg->setColor(c,i, j);

        }
    }

    outimg->flipVertical();

    if(!outimg) return 0;
    if(needResample){
        osg::Image* img2 = new osg::Image();
        img2->allocateImage(img_size,img_size,4,GL_RGBA,GL_UNSIGNED_BYTE);
        mapImage(outimg,ext2,img2,extent);
        return img2;
    }else{
        return outimg;
    }
}

// From easyrgb.com
float Hue_2_RGB( float v1, float v2, float vH )
{
    if ( vH < 0 ) vH += 1;
    if ( vH > 1 ) vH -= 1;
    if ( ( 6 * vH ) < 1 ) return ( v1 + ( v2 - v1 ) * 6 * vH );
    if ( ( 2 * vH ) < 1 ) return ( v2 );
    if ( ( 3 * vH ) < 2 ) return ( v1 + ( v2 - v1 ) * ( ( 2 / 3 ) - vH ) * 6 );
    return ( v1 );
}


osg::Image *DataSourceImageGDAL::__read(int windowX, int windowY, int windowWidth, int windowHeight, int destWidth, int destHeight)
{

    int numBands = poDataset->GetRasterCount();


    GDALRasterBand* bandGray = 0;
    GDALRasterBand* bandRed = 0;
    GDALRasterBand* bandGreen = 0;
    GDALRasterBand* bandBlue = 0;
    GDALRasterBand* bandAlpha = 0;
    GDALRasterBand* bandPalette = 0;

    int internalFormat = GL_LUMINANCE;
    unsigned int pixelFormat = GL_LUMINANCE;
    unsigned int dataType = 0;
    unsigned int numBytesPerComponent = 0;

    GDALDataType targetGDALType = GDT_Byte;

    for(int b=1;b<=numBands;++b)
    {

        GDALRasterBand* band = poDataset->GetRasterBand(b);
        if (band->GetColorInterpretation()==GCI_GrayIndex) bandGray = band;
        else if (band->GetColorInterpretation()==GCI_RedBand) bandRed = band;
        else if (band->GetColorInterpretation()==GCI_GreenBand) bandGreen = band;
        else if (band->GetColorInterpretation()==GCI_BlueBand) bandBlue = band;
        else if (band->GetColorInterpretation()==GCI_AlphaBand) bandAlpha = band;
        else if (band->GetColorInterpretation()==GCI_PaletteIndex) bandPalette = band;
        else bandGray = band;

        if (dataType==0)
        {
            targetGDALType = band->GetRasterDataType();
            switch(band->GetRasterDataType())
            {
            case(GDT_Byte): dataType = GL_UNSIGNED_BYTE; numBytesPerComponent = 1; break;
            case(GDT_UInt16): dataType = GL_UNSIGNED_SHORT; numBytesPerComponent = 2; break;
            case(GDT_Int16): dataType = GL_SHORT; numBytesPerComponent = 2; break;
            case(GDT_UInt32): dataType = GL_UNSIGNED_INT; numBytesPerComponent = 4; break;
            case(GDT_Int32): dataType = GL_INT; numBytesPerComponent = 4; break;
            case(GDT_Float32): dataType = GL_FLOAT; numBytesPerComponent = 4; break;
            case(GDT_Float64): dataType = GL_DOUBLE; numBytesPerComponent = 8; break;  // not handled
            default: dataType = 0; numBytesPerComponent = 0; break; // not handled
            }
        }
    }

    int s = destWidth;
    int t = destHeight;
    int r = 1;

    if (dataType==0)
    {
        dataType = GL_UNSIGNED_BYTE;
        numBytesPerComponent = 1;
        targetGDALType = GDT_Byte;
    }

    unsigned char* imageData = 0;

    if (bandRed && bandGreen && bandBlue)
    {
        if (bandAlpha)
        {
            // RGBA

            int pixelSpace=4*numBytesPerComponent;
            int lineSpace=destWidth * pixelSpace;

            imageData = new unsigned char[destWidth * destHeight * pixelSpace];
            pixelFormat = GL_RGBA;
            internalFormat = GL_RGBA;

            OSG_INFO << "reading RGBA"<<std::endl;

            bandRed->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandGreen->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+1*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandBlue->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+2*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandAlpha->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+3*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);

        }
        else
        {
            // RGB

            int pixelSpace=3*numBytesPerComponent;
            int lineSpace=destWidth * pixelSpace;

            imageData = new unsigned char[destWidth * destHeight * pixelSpace];
            pixelFormat = GL_RGB;
            internalFormat = GL_RGB;

            OSG_INFO << "reading RGB"<<std::endl;

            bandRed->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandGreen->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+1*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandBlue->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+2*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);

        }
    }
    else if (bandGray)
    {
        if (bandAlpha)
        {
            // Luminance alpha
            int pixelSpace=2*numBytesPerComponent;
            int lineSpace=destWidth * pixelSpace;

            imageData = new unsigned char[destWidth * destHeight * pixelSpace];
            pixelFormat = GL_LUMINANCE_ALPHA;
            internalFormat = GL_LUMINANCE_ALPHA;

            OSG_INFO << "reading grey + alpha"<<std::endl;

            bandGray->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandAlpha->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+1*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
        }
        else
        {
            // Luminance map
            int pixelSpace=1*numBytesPerComponent;
            int lineSpace=destWidth * pixelSpace;

            imageData = new unsigned char[destWidth * destHeight * pixelSpace];
            pixelFormat = GL_LUMINANCE;
            internalFormat = GL_LUMINANCE;

            OSG_INFO << "reading grey"<<std::endl;

            bandGray->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
        }
    }
    else if (bandAlpha)
    {
        // alpha map
        int pixelSpace=1*numBytesPerComponent;
        int lineSpace=destWidth * pixelSpace;

        imageData = new unsigned char[destWidth * destHeight * pixelSpace];
        pixelFormat = GL_ALPHA;
        internalFormat = GL_ALPHA;

        OSG_INFO << "reading alpha"<<std::endl;

        bandAlpha->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);

    }
    else if (bandPalette)
    {
        // Paletted map
        int pixelSpace=1*numBytesPerComponent;
        int lineSpace=destWidth * pixelSpace;

        unsigned char *rawImageData;
        rawImageData = new unsigned char[destWidth * destHeight * pixelSpace];
        imageData = new unsigned char[destWidth * destHeight * 4/*RGBA*/];
        pixelFormat = GL_RGBA;
        internalFormat = GL_RGBA;

        OSG_INFO << "reading palette"<<std::endl;
        OSG_INFO << "numBytesPerComponent: " << numBytesPerComponent << std::endl;

        bandPalette->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(rawImageData),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);

        // Map the indexes to an actual RGBA Value.
        for (int i = 0; i < destWidth * destHeight; i++)
        {
            const GDALColorEntry *colorEntry = bandPalette->GetColorTable()->GetColorEntry(rawImageData[i]);
            GDALPaletteInterp interp = bandPalette->GetColorTable()->GetPaletteInterpretation();
            if (!colorEntry)
            {
                //FIXME: What to do here?

                //OSG_INFO << "NO COLOR ENTRY FOR COLOR " << rawImageData[i] << std::endl;
                imageData[4*i+0] = 255;
                imageData[4*i+1] = 0;
                imageData[4*i+2] = 0;
                imageData[4*i+3] = 1;

            }
            else
            {
                if (interp == GPI_RGB)
                {
                    imageData[4*i+0] = colorEntry->c1;
                    imageData[4*i+1] = colorEntry->c2;
                    imageData[4*i+2] = colorEntry->c3;
                    imageData[4*i+3] = colorEntry->c4;
                }
                else if (interp == GPI_CMYK)
                {
                    // from wikipedia.org
                    short C = colorEntry->c1;
                    short M = colorEntry->c2;
                    short Y = colorEntry->c3;
                    short K = colorEntry->c4;
                    imageData[4*i+0] = 255 - C*(255 - K) - K;
                    imageData[4*i+1] = 255 - M*(255 - K) - K;
                    imageData[4*i+2] = 255 - Y*(255 - K) - K;
                    imageData[4*i+3] = 255;
                }
                else if (interp == GPI_HLS)
                {
                    // from easyrgb.com
                    float H = colorEntry->c1;
                    float S = colorEntry->c3;
                    float L = colorEntry->c2;
                    float R, G, B;
                    if ( S == 0 )                       //HSL values = 0 - 1
                    {
                        R = L;                      //RGB results = 0 - 1
                        G = L;
                        B = L;
                    }
                    else
                    {
                        float var_2, var_1;
                        if ( L < 0.5 )
                            var_2 = L * ( 1 + S );
                        else
                            var_2 = ( L + S ) - ( S * L );

                        var_1 = 2 * L - var_2;

                        R = Hue_2_RGB( var_1, var_2, H + ( 1 / 3 ) );
                        G = Hue_2_RGB( var_1, var_2, H );
                        B = Hue_2_RGB( var_1, var_2, H - ( 1 / 3 ) );
                    }
                    imageData[4*i+0] = static_cast<unsigned char>(R*255.0f);
                    imageData[4*i+1] = static_cast<unsigned char>(G*255.0f);
                    imageData[4*i+2] = static_cast<unsigned char>(B*255.0f);
                    imageData[4*i+3] = static_cast<unsigned char>(255.0f);
                }
                else if (interp == GPI_Gray)
                {
                    imageData[4*i+0] = static_cast<unsigned char>(colorEntry->c1*255.0f);
                    imageData[4*i+1] = static_cast<unsigned char>(colorEntry->c1*255.0f);
                    imageData[4*i+2] = static_cast<unsigned char>(colorEntry->c1*255.0f);
                    imageData[4*i+3] = static_cast<unsigned char>(255.0f);
                }
            }
        }
        delete [] rawImageData;
    }
    else
    {
        OSG_INFO << "not found any usable bands in file."<<std::endl;
    }


    if (imageData)
    {
        osg::Image* image = new osg::Image;
        //image->setFileName(fileName.c_str());
        image->setImage(s,t,r,
                        internalFormat,
                        pixelFormat,
                        dataType,
                        (unsigned char *)imageData,
                        osg::Image::USE_NEW_DELETE);

        image->flipVertical();

        return image;
    }
    return 0;
}

osg::Image *DataSourceImageGDAL::__reads(int windowX, int windowY, int windowWidth, int windowHeight, int destWidth, int destHeight)
{
    int numBands = poDataset->GetRasterCount();


    GDALRasterBand* bandGray = 0;
    GDALRasterBand* bandRed = 0;
    GDALRasterBand* bandGreen = 0;
    GDALRasterBand* bandBlue = 0;
    GDALRasterBand* bandAlpha = 0;
    GDALRasterBand* bandPalette = 0;

    int internalFormat = GL_LUMINANCE;
    unsigned int pixelFormat = GL_LUMINANCE;
    unsigned int dataType = 0;
    unsigned int numBytesPerComponent = 0;

    GDALDataType targetGDALType = GDT_Byte;

    for(int b=1;b<=numBands;++b)
    {

        GDALRasterBand* band = poDataset->GetRasterBand(b);
        if (band->GetColorInterpretation()==GCI_GrayIndex) bandGray = band;
        else if (band->GetColorInterpretation()==GCI_RedBand) bandRed = band;
        else if (band->GetColorInterpretation()==GCI_GreenBand) bandGreen = band;
        else if (band->GetColorInterpretation()==GCI_BlueBand) bandBlue = band;
        else if (band->GetColorInterpretation()==GCI_AlphaBand) bandAlpha = band;
        else if (band->GetColorInterpretation()==GCI_PaletteIndex) bandPalette = band;
        else bandGray = band;

        if (dataType==0)
        {
            targetGDALType = band->GetRasterDataType();
            switch(band->GetRasterDataType())
            {
            case(GDT_Byte): dataType = GL_UNSIGNED_BYTE; numBytesPerComponent = 1; break;
            case(GDT_UInt16): dataType = GL_UNSIGNED_SHORT; numBytesPerComponent = 2; break;
            case(GDT_Int16): dataType = GL_SHORT; numBytesPerComponent = 2; break;
            case(GDT_UInt32): dataType = GL_UNSIGNED_INT; numBytesPerComponent = 4; break;
            case(GDT_Int32): dataType = GL_INT; numBytesPerComponent = 4; break;
            case(GDT_Float32): dataType = GL_FLOAT; numBytesPerComponent = 4; break;
            case(GDT_Float64): dataType = GL_DOUBLE; numBytesPerComponent = 8; break;  // not handled
            default: dataType = 0; numBytesPerComponent = 0; break; // not handled
            }
        }
    }

    int s = destWidth;
    int t = destHeight;
    int r = 1;

    if (dataType==0)
    {
        dataType = GL_UNSIGNED_BYTE;
        numBytesPerComponent = 1;
        targetGDALType = GDT_Byte;
    }

    unsigned char* imageData = 0;

    if (bandRed && bandGreen && bandBlue)
    {
        if (bandAlpha)
        {
            // RGBA

            int pixelSpace=4*numBytesPerComponent;
            int lineSpace=destWidth * pixelSpace;

            imageData = new unsigned char[destWidth * destHeight * pixelSpace];
            pixelFormat = GL_RGBA;
            internalFormat = GL_RGBA;

            OSG_INFO << "reading RGBA"<<std::endl;

            bandRed->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandGreen->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+1*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandBlue->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+2*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandAlpha->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+3*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);

        }
        else
        {
            // RGB

            int pixelSpace=3*numBytesPerComponent;
            int lineSpace=destWidth * pixelSpace;

            imageData = new unsigned char[destWidth * destHeight * pixelSpace];
            pixelFormat = GL_RGB;
            internalFormat = GL_RGB;

            OSG_INFO << "reading RGB"<<std::endl;

            bandRed->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandGreen->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+1*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandBlue->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+2*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);

        }
    }
    else if (bandGray)
    {
        if (bandAlpha)
        {
            // Luminance alpha
            int pixelSpace=2*numBytesPerComponent;
            int lineSpace=destWidth * pixelSpace;

            imageData = new unsigned char[destWidth * destHeight * pixelSpace];
            pixelFormat = GL_LUMINANCE_ALPHA;
            internalFormat = GL_LUMINANCE_ALPHA;

            OSG_INFO << "reading grey + alpha"<<std::endl;

            bandGray->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
            bandAlpha->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+1*numBytesPerComponent),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
        }
        else
        {
            // Luminance map
            int pixelSpace=1*numBytesPerComponent;
            int lineSpace=destWidth * pixelSpace;

            imageData = new unsigned char[destWidth * destHeight * pixelSpace];
            pixelFormat = GL_LUMINANCE;
            internalFormat = GL_LUMINANCE;

            OSG_INFO << "reading grey"<<std::endl;

            bandGray->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);
        }
    }
    else if (bandAlpha)
    {
        // alpha map
        int pixelSpace=1*numBytesPerComponent;
        int lineSpace=destWidth * pixelSpace;

        imageData = new unsigned char[destWidth * destHeight * pixelSpace];
        pixelFormat = GL_ALPHA;
        internalFormat = GL_ALPHA;

        OSG_INFO << "reading alpha"<<std::endl;

        bandAlpha->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(imageData+0),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);

    }
    else if (bandPalette)
    {
        // Paletted map
        int pixelSpace=1*numBytesPerComponent;
        int lineSpace=destWidth * pixelSpace;

        unsigned char *rawImageData;
        rawImageData = new unsigned char[destWidth * destHeight * pixelSpace];
        imageData = new unsigned char[destWidth * destHeight * 4/*RGBA*/];
        pixelFormat = GL_RGBA;
        internalFormat = GL_RGBA;

        OSG_INFO << "reading palette"<<std::endl;
        OSG_INFO << "numBytesPerComponent: " << numBytesPerComponent << std::endl;

        bandPalette->RasterIO(GF_Read,windowX,windowY,windowWidth,windowHeight,(void*)(rawImageData),destWidth,destHeight,targetGDALType,pixelSpace,lineSpace);

        // Map the indexes to an actual RGBA Value.
        for (int i = 0; i < destWidth * destHeight; i++)
        {
            const GDALColorEntry *colorEntry = bandPalette->GetColorTable()->GetColorEntry(rawImageData[i]);
            GDALPaletteInterp interp = bandPalette->GetColorTable()->GetPaletteInterpretation();
            if (!colorEntry)
            {
                //FIXME: What to do here?

                //OSG_INFO << "NO COLOR ENTRY FOR COLOR " << rawImageData[i] << std::endl;
                imageData[4*i+0] = 255;
                imageData[4*i+1] = 0;
                imageData[4*i+2] = 0;
                imageData[4*i+3] = 1;

            }
            else
            {
                if (interp == GPI_RGB)
                {
                    imageData[4*i+0] = colorEntry->c1;
                    imageData[4*i+1] = colorEntry->c2;
                    imageData[4*i+2] = colorEntry->c3;
                    imageData[4*i+3] = colorEntry->c4;
                }
                else if (interp == GPI_CMYK)
                {
                    // from wikipedia.org
                    short C = colorEntry->c1;
                    short M = colorEntry->c2;
                    short Y = colorEntry->c3;
                    short K = colorEntry->c4;
                    imageData[4*i+0] = 255 - C*(255 - K) - K;
                    imageData[4*i+1] = 255 - M*(255 - K) - K;
                    imageData[4*i+2] = 255 - Y*(255 - K) - K;
                    imageData[4*i+3] = 255;
                }
                else if (interp == GPI_HLS)
                {
                    // from easyrgb.com
                    float H = colorEntry->c1;
                    float S = colorEntry->c3;
                    float L = colorEntry->c2;
                    float R, G, B;
                    if ( S == 0 )                       //HSL values = 0 - 1
                    {
                        R = L;                      //RGB results = 0 - 1
                        G = L;
                        B = L;
                    }
                    else
                    {
                        float var_2, var_1;
                        if ( L < 0.5 )
                            var_2 = L * ( 1 + S );
                        else
                            var_2 = ( L + S ) - ( S * L );

                        var_1 = 2 * L - var_2;

                        R = Hue_2_RGB( var_1, var_2, H + ( 1 / 3 ) );
                        G = Hue_2_RGB( var_1, var_2, H );
                        B = Hue_2_RGB( var_1, var_2, H - ( 1 / 3 ) );
                    }
                    imageData[4*i+0] = static_cast<unsigned char>(R*255.0f);
                    imageData[4*i+1] = static_cast<unsigned char>(G*255.0f);
                    imageData[4*i+2] = static_cast<unsigned char>(B*255.0f);
                    imageData[4*i+3] = static_cast<unsigned char>(255.0f);
                }
                else if (interp == GPI_Gray)
                {
                    imageData[4*i+0] = static_cast<unsigned char>(colorEntry->c1*255.0f);
                    imageData[4*i+1] = static_cast<unsigned char>(colorEntry->c1*255.0f);
                    imageData[4*i+2] = static_cast<unsigned char>(colorEntry->c1*255.0f);
                    imageData[4*i+3] = static_cast<unsigned char>(255.0f);
                }
            }
        }
        delete [] rawImageData;
    }
    else
    {
        OSG_INFO << "not found any usable bands in file."<<std::endl;
    }


    if (imageData)
    {
        osg::Image* image = new osg::Image;
        //image->setFileName(fileName.c_str());
        image->setImage(s,t,r,
                        internalFormat,
                        pixelFormat,
                        dataType,
                        (unsigned char *)imageData,
                        osg::Image::USE_NEW_DELETE);

        //        image->flipVertical();

        return image;
    }
    return 0;
}


